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Abstract 

With the imminent start of LHC experiments, development of phenomenological 
tools, and in particular the Monte Carlo programs and algorithms, becomes urgent. 
A new algorithm for the generation of a parton shower initiated by the single ini- 
tial hadron beam is presented. The new algorithm is of the class of the so called 
"constrained MC" type algorithm (an alternative to the backward evolution MC 
algorithm) , in which the energy and the type of the parton at the end of the parton 
shower are constrained (predefined). The complete kinematics configurations with 
explicitly constructed four momenta are generated and tested. Evolution time is 
identical with rapidity and minimum transverse momentum is used as an infrared 
cut-off. All terms of the leading-logarithmic approximation in the DGLAP evolu- 
tion are properly accounted for. In addition, the essential improvements towards 
the so-called CCFM/BFKL models are also properly implemented. The resulting 
parton distributions are cross-checked up to the 10 -3 precision level with the help 
of a multitude of comparisons with other MC and non-MC programs. We regard 
these tests as an important asset to be exploited at the time when the presented MC 
will enter as a building block in a larger MC program for W/Z production process 
at LHC. 
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1 Introduction 



In the past, as at present, the central goal of high energy physics is to explore new ranges 
of energies of interactions. These new ranges of energies either facilitate the discovery of 
new particles and interactions or validate the Standard Model of elementary interactions, 
as understood today, by extending them to an even broader range of energies (distances) 
than currently available. The new generation of experiments at the nearly completed 
Large Hadron Collider (LHC) in Geneva will be ready soon to take data. 

For the proper interpretation of expected new data, an effort in understanding known 
physics is needed. In particular, it might be that the signatures of the new physics will 
have to be deciphered from the background of dominant processes expected from the 
Standard Model. 

In the case of hadron colliders, description of the Standard Model processes is rather 
complicated; the colliding hadron beams are not the elementary fields of the Standard 
Model but bounded states of quarks and gluons. Even worse, in the low energy limit quan- 
tum Chromodynamics (QCD) looses predictive power and does not control the relations 
between the hadron wave function and elementary fields representing quarks and gluons. 
It is necessary, albeit highly nontrivial, to combine a phenomenological description of low 
energy strong interaction phenomena with the rigorous perturbative QCD predictions at 
high energies. 

A multitude of techniques have been developed to merge low energy aspects of strong 
interaction with the high energy calculations from perturbative QCD. In this work we con- 
centrate on the methodology based on the so-called parton distribution functions (PDF) 
and parton shower Monte Carlo (PSMC). Special attention will be payed to technical 
aspects, in particular to precision testing of the numerical tools. We believe that this is 
very important for the future efforts in minimizing overall systematic errors of the QCD 
prediction in which PSMCs are used. 

In this work we shall concentrate on the question of the evolution equation of the PDF 
in the Monte Carlo (MC) form. Such an evolution equation for the initial state hadron 
describes how the PDF responds to an increase of the dimensional, large energy scale 
Q = /i, set by the hard process probing the PDF. The formula for the integrated cross 
section, which combines the hard process with the matrix element, has been proved within 
perturbative QCD in a form of the so called factorization theorems, see for instance refs. 
[TJ [2] . These theorems have been proven starting directly from the Feynman diagrams, 
integrated over the phase space and convoluted with the nonperturbative parton wave 
function in a hadron. The evolution equation of the PDF can also be formulated using 
renormalization group and operator product expansion [3]. 

However, for the real-life practice of the present and future hardon-hadron and hadron- 
electron collider experiments, one needs a more refined (exclusive) picture of the multi- 
parton production, the simplest one being the so called parton shower, governed by per- 
turbative QCD. In principle, it should reproduce the evolution of PDFs, after integrating 
it over the Lorentz invariant phase space. If the above is conveniently implemented in 
the form of the parton shower Monte Carlo event generator, see for example [HE], with 
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the inclusion of the parton hadronization process - a very useful feature for the collider 
experiments. In the construction of such a PSMC the evolution equation of the PDF is 
used as a guide for defining distributions of the partons emitted from a single energetic 
hadron (the initiatial parton in the shower )0. 

Having in mind the above context, this work has several aims. The principal aim is to 
use once again the evolution equation of the PDF in order to model the multiparton parton 
shower initiated by the single parton located inside the single hadron of the collider beam. 
We shall insist that, as in the factorization theorems, this modelling has to be universal, 
that is independent from showering of the other hadron beam, spectator partons, and the 
type of the hard process of the parton-parton scattering at the large energy scale Q. 

Another aim is to define and maintain a clear prescription relating variables in the 
evolution of the PDF with the four-momenta in the PSMC. We keep in mind that such 
a PSMC will be a building block to be used for two beams^l, hence the upper limit of the 
multiparton phase space (related to Q) should allow for smooth coverage of the entire 
phase space, without any gaps and overlaps. 

There is also an important technical problem to be addressed: the off-shell parton 
entering the hard process has to have predefined energy and flavour matching preferences 
of the hard process; hence, a Markovian MC (which is a natural MC implementation of the 
PSMC) cannot be used to model the initial state PSMC. However, instead of using the so 
called backward evolution [7], our choice will be to employ the technique of the constrained 
MC (referred to as CMC technique); that is to generate the multiparton distribution with 
the restriction on the value of the parton energy and type of the parton. Two distinct 
versions of this relatively new technique were proposed and tested in refs. [El El [10] for 
the DGLAP |TTJ evolution in the leading-logarithmic (LL) approximation. 

The aim of this work is to extend the most promising variant of the above CMC 
technique to a wider class of evolution kernels, beyond DGLAP, towards evolution models 
of the CCFM class [12]; maintaining at the same time the explicit mapping of the evolution 
variables into four-momenta. 

The other important longer term goal is to facilitate the inclusion of the complete 
NLO corrections by means of the clearer /cleaner modelling of the PSMC, as compared 
with the existing combined NLO and PSMC calculations [131 [H]. This will be achieved, 
for example, by means of better coverage of the phase space in the basic parton shower 
MC. 

The outline of the paper is the following: In section 2 we shall formulate the general 
formalism of the evolution equations and their solutions in a form suitable for the CMC 
technique. In section 3 we discuss two particular types of the kernels and related Sudakov 
form-factors. In section 4 we outline the CMC algorithm for the pure gluonstrahlung 
segments. In section 5 we shall introduce in the CMC quark gluon transitions. In section 
6 results of precision numerical tests of our CMC implementation will be reported. For 

lr This is clearly a kind of "backward engineering" - it would be better to get distributions of partons 
forming PDF at large Q directly from the Feyman diagrams. Unfortunately it is too difficult. 

2 This will be done in the forthcomming paper on the new parton shower MC for W/Z production in 
hadron collider of the forthcomming paper [6]. 
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testing CMC implementations we shall use auxiliary Markovian MC programs which are 
described and tested in separate papers [15] . Finally, we summarize the main results; 
some technical details will be included in the appendices. 



2 Evolution equations and solutions 

In the following we shall formulate the mathematical framework for the evolution equation 
and its solutions in a form suitable for the construction of the CMC algorithm in the latter 
part of the paper. 

The generic evolution equation covering several types of evolution reads 

d t D f (t,x) = [ duX ff (t,x,u)D f (t,u). (1) 

fl J x 

It describes the evolution of the parton distribution function Dj(t, u), where x is fraction 
of the hadron momentum carried by the parton and j is the type (flavour) of the parton. 
The variable t = In Q is traditionally called an evolution time and it represents the (large) 
energy scale Q — \i at which the PDF is probed using a hard scattering process. The LL 
DGLAP case [11] is covered by eq. ([T]) with the following identification 

x\ = a s (t) 2 p(0 ) / or 
v uJ 2% u J* \ ' u> 



X ff , (t,x,u) = - OV ( t, - ) = - P% ( t, - ) , (2) 



where P^,(z) is the standard LL DGLAP kernel and the factor 2 is related to our definition 
of the evolution variable t. 

In the compact operator (matrix) notation eq. ([1]) reads 

dtD(t) = K(t) D(t). (3) 

Given a known D(to), the formal solution at any later "time" t > to is provided by the 
time ordered exponential 

D(t) = exp (y K(t')dt'^ D(t ) = G K (t,t ) D(t ). (4) 
The time- ordered exponential evolution operator readg^l 

/ K(t')dt'\ =l + J2\[ <ttA>U-iK(ti), (5) 

J to J T n= i i= i J t Q 

For the sake of completeness, let us write the explicit definition^ of the multiplication 
operation as used and defined in eqs. ([3][5j: 

(K(t 2 ) K(ti)) /2ji (x 2 ,xi) = Y dx ' % hf'( t ^ x ^ x ') Xf'h(ti,a?,xi). (6) 

fl J X2 



3 Here and in the following we define Il"=i A* = -^n-^n-i ■ ■ ■ A2A1. 

4 In the case of QCD evolution K(t,) transforms u G (0,1) into x obeying x < u. This is due to 
4-momentum conservation. 
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We shall often be dealing with the case of the kernel split into parts, for example: 

K(t) =K A (t) + K B (t). (7) 
In such a case the solution of eq. (jl]) can be reorganized as follow^ 



D(i) = G K B(Mo)D(io) + ^ 

71=1 
OO 

= G K B(t,to)D(t ) + X) 



[ dti 9 ti>ti l G K B(t i+1 , ti) K A (ti 



71=1 



1=1 " *o 



G KS (^,to)D(t ) 



' / dU 

i=i J ti-i 



G K s(t, t n ) 



jK A (tj) G K s(tj,tj_x^ 



i=l 



D(t ) 



where t n+ i = £ and Grb is the evolution operator of eq. (j5J) of the evolution with the 
kernel ~K. B . Formal proof of eq. (jHJ) is given in ref. [T6] . 



<n-l4-l 



n+1 11 j iA j 

t. . f. 



1111-1 i r l o 1 



*n-l 



Figure 1: The scheme of integration variables and summation indices in eq. ( flOl) : the circles 
correspond to K A (tj) and the ovals to G k b, z[ = x'Jxi-i and Zi = Xijx\. 

In the following we are going to nest eq. (jSJ) twice. First, we employ it in order to 
isolate gluonstrahlung and flavour- changing parts of the evolution, exploiting the following 
split of the kernel 

X(t) ff , = X A (t) ff + X B (t) ff , = (1 - 5 ff )X(t) ff + 5 ff X(t) ff . (9) 
In this case G K s represents pure gluonstrahlung and is diagonal in the flavour index. In 



5 Here and in the following, we understand that the scope of the indices ceases at the closing bracket. 
However, the validity scope of indiced variables, for instance of ti, extends until the formula's end. The 
use of eq. ((6]) is understood to be adjusted accordingly. 
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the standard integro-tensorial notation eq. (jHJ) looks as follows: 

\ oo 

D f (t,x)= dx G ff (K B ] t, to; x, x ) D f (t , x ) + 

rp n=l / n _i,. 



/o 



X 



X 



r 

y n y ^ y ^ y ^ ^^^-i 

~ 2 = 1 + ' n n 



ti_i o 



G ff (K B ;t,t n ,x,x n ) (10) 



i=l 



Df (t ,X ) 



where / n = /. In the above and the following equations we adopt the following notation: 



5(x - y) 



and 



0*, 



/ j/<x — 1 for y < x and 9 y<x = for y > x. 

Similarly, 8 z<y<x = 6 z<y 6 y<x . The chain of integration variables and flavour indices is 
depicted schematically in Fig. [TJ 

Next, eq. (jHJ) is used in order to resum the virtual IR-divergent part % v of the gluon- 
strahlung kernel % B 



X B f ,(t,x,u) = 5ff'X f f(t,x,u) = 6f f > (pCf f {t,x, u) + Xf f (t,x, u)), 
X^j{t, x, u) = —5 X=U X j jr(t, x), 

Xff(t,X,u) = 9 x <u-A(x,u)X^ff(t 



X, u 



where A(x,u) is finite IR cut-off, not necessarily infinitesimal. In order to resum (expo- 
nentiate) the virtual part X v of the kernel the following version of eq. (jHJ) is employed 



Gkv+k«(Mo) = G K v-(t, to) + 



71=1 



n pt 

' / dt i e ti>t ._ 1 G Ji v(t i+1 ,t i )K R (t i ) 
i=l 



G K v(ti, to)- 



(12) 



Since K (t«) is diagonal in x, u and in the flavour index, and also because of 



{G K v{t i+1 ,ti)} ff {x,u) = S x=u g-*/^!.**!*), $^+1,^) = 

we obtain immediately 

G ff (K B ]t b ,t a ,x,u) = {G K B(t b ,t a )}ff(x,u) = 



dtX v ff (t 1 x), (13) 



n=l L i=l " ta 
x e -*/(*6>*n|a;) 



(14) 



Jac? / (t iI x il x < _i) e - ft /*- t '- i i-«- i J 



where t n+ i = to = t a and xq = u. The above result can also be obtained by iterating 
the evolution equation for G K s(t,to) with the boundary condition G k b(£o,£o) = I, see 
for instance ref. |17j . 

The algebra resulting from (jSJ) is quite general and does not rely on any particular 
form of the kernel. For example, in the above calculations we did not have to invoke the 
energy sum rules, or any other specific restrictions on the overall normalization embodied 
usually in the virtual corrections. Also, we did not define yet the relation between the 
evolution variables U and Xj and parton four-momenta. This will be done in the following 
section. 

3 Evolution kernel and variables 

Before we specify details of the evolution kernels used in this work, let us discuss the 
relation between the evolution variables Xj, tj and the emitted parton four-momenta kj 
in the corresponding parton shower MC. 

In the proofs of the factorization theorems [TJ [18j [2] and its practical realizations 
like in ref. [TH], typically within MS scheme, one projects the four-momenta of the (off- 
shell) partons into the 1-dimensional variable of the evolution, typically the dimensionless 
lightcone variable x. The so-called factorization scale fip = Q measuring the size of the 
available parton emission phase space is usually set by the kinematics of the hard pro- 
cess. The underlying QCD differential distribution (the QCD matrix element times the 
phase space) is reduced to a chain of parton splittings with the x variable of the PDF 
being the fraction of the initial hadron energy Eh carried by the parton entering the hard 
process. The variable e* = Q defines the boundary (maximum value) for many ordered 
variables^) U. The variable t may be related directly to one of the phase space parame- 
ters (virtuality, transverse momentum, angle) or the abstract dimensional scale variable 
\ip resulting from the formal procedure of cancelling IR singularities in the dimensional 
regularization method. In the classical construction of the parton shower MC, one must 
invert the above mapping of the phase space variables into the evolution variables; that 
is to construct parton four-momenta out of ti and Xi and to reconstruct fully differential 
parton distributions in terms of these four-momenta. Obviously this procedure is not 
unique and requires some guidance from the detailed knowledge of the structure of the 
IR singularities 3 of the original QCD matrix element. 

In this work we shall construct the CMC algorithm for two new types of generalized 
kernels, in addition to the ones of LL DGLAP of eq. (j2J) for which examples of CMC were 
already constructed in refs. [8] and [9]. The new kernels are based on the following generic 
form 




K. ffl (t,x,u) = K.ff(t,x,u)6 u > x+A (t,u), 



(15) 



6 This simplified picture is valid at least in the leading-logarithmic approximation. 
7 We understand IR singularities as both the collinear and soft ones. 
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Figure 2: Two types of the infrared (IR) boundaries on the u — x plane (a) u — x > ue and 
(d) u — x > e for the triangle (red) which depicts the area where the real emission part K, e is 
nonzero. The diagonal line (blue) represents places where the virtual K? is nonzero. 

where P^,(t, z) is the standard LL kernel (DGLAP) and t-dependence enters into it only 
through the IR regulator A(t,u). The new kernel types correspond to different choices 
of the argument of the strong coupling constant and to different forms of the A(t,u) 
regulator. The two new types of the regulator A(t, u) used in this work are depicted on 
the x — u plane in Fig. [2] and will be defined in the following section. 

The other important departure from DGLAP is the strong coupling constant cts(t, x, u) 
which now may depend on all evolution variables. We shall consider the strong coupling 
constant a$ depending on z = x/u or on the transverse momentum k T defined below. 
Before we define the evolution kernel in a detail, we have to elaborate first the mapping 
of the evolution variables into four-momenta. 

3.1 Relating evolution variables to four- moment a 

The essential decision in the construction of the parton shower MC concerns the choice of 
a kinematics variable in the solutions of the evolution equations ffTOl) and (IT%|) ; this choice 
is in one-to-one correspondence with the evolution time variable ti and the limiting value 
t. We choose to associate t with the rapidity (angle) of the emitted parton, following the 
well known arguments on the colour coherence exposed in many papers, see for instance 
refs. [201 Ell [22] and further references therein. 

We define the lightcone variables = q° ± q 3 and normalize the parton momenta 
with respect to the energy E h of the initial massless hadron, see Fig. [31 

q£ = 2Eh, and qf = Xi2Eh, where < Xi < 1. (16) 

These relations hold in the rest frame of the hard process system (HRS) with the z-axis 
along the momentum of initial state hadron h. 
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Figure 3: Kinematics in the evolution history (tree). 

In particular, for the parton initiating a parton cascade, we have 

q+ = x 2E h . (17) 

This parton has negligible transverse momentunjfl, k$ — 0. In the HRS each on-shell i-th 
emitted particle will take away a part of the lightcone variable 

H = qt-i - vt = (x>-i - Xi)2E h . (is) 

The above is not enough to define k^. For this we need to define at least k T = \k T \ or 
the rapidity 77$ = \ hx{kf /k^). Once the azimuthal angle ipi is added, we can complete 
the mapping from the evolution variables to the 4-momentum of the i-th emitted parton 

kj^ {tii %ii 1) ^Pi)- 

If we associate the evolution time ti with the rapidity rji then the following relation 




KK = K\hT = fc i e ~* = (Zi-i -x t )2E h e~^\ (19) 



valid for kf = 0, provides the transverse momentum and thus fcf . We can also eliminate 
kj with the help of the following conventional relation which defines the evolution time t 

kf Ee^n-i-ii). (20) 

Note that this relation translates into k T = e l (u — x) in the general evolution equation 
of eq. ([T]). We observe that evolution time ti and rapidity 77, are related by a linear 
transformation of the following explicit form: 

e ti = e H2E h )- m ^ Vi = \ n (2E h ) - ti. (21) 

The above relation is the main result of this section. 

Summarizing, the mapping of ti and Xi into 4-momenta kf,i = 1,2, ...n, can be 
written now in an explicit manner: 



K = (xi-i - Xi)2E h , kf = - Xi )e u , k, = {kf y/k?, 

k° = (kt + K)/2, k* = (kt-K)/2. [ } 



3 In the realistic MC k$ will be distributed according to a Gaussian profile with the width of ~ 0(X). 
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Last but not least, we have to define also the phase space limits, U G (t m m,t max ) 
and rji G (r) m i a , r/ max )- One has to be very careful at this step. The maximum evolution 
time t max (minimum rapidity T] m[n ) is set by the requirement that all emitted partons are 
confined to the forward hemispher^], 90° > 6*j, or equivalently rji > ?7 m i n = 0. This implies 

t max = ln(2£ h ), E h = ie w . (23) 

The minimum evolution time is determined by the phase space opening point for the 
first emission, due to minimum transverse momentum kf in = A: 

e h x > A. (24) 

This leads to 

ti> t x — lnx , t\ = In A 

and therefore to 

t>t n > t n _! > . . . > t 2 > h > t x - In x = t min . 
This automatically determines the maximum rapidity 

?7max = ln(2J5 h ) - In A + lnx . 

Altogether 

U G (tmm,imax) = (In A - In x Q , ]n(2E h )) , 

Vi e (^min,^max) = (0, ln(2E h ) + lnx - In A)), (25) 

^7max ^7min ^max ^min ln(2£'/ l / A). 

Let us remark that the naive assignment t max = ln(^), without the factor of 2, would 
lead, because of eq. (12~TI) . to a partial coverage of the forward hemisphere only, rji > ln(2). 
On the other hand, this factor of 2 may look justified; the absolute kinematic range of 
the transverse momentum is kf G (A, fc^), where fc^ ax = e* max = 2^ results from the 
relation kf = e**(xj_i — Xj) and energy conservation Xi ^ 1. The reader may notice that 
this limit is a factor of 2 higher than in the familiar inequality kf < resulting from 
the 4-momentum conservation operating in both hemispheres simultaneously, i.e. on kf 
and k~. However, our limit is valid, including the factor of 2, for a single hemisphere 
separated from any "activity" on the other side! 

Finally, we illustrate the real emission phase space in Fig. HJ[a) using the rapidity 
variable r\ and the log of transverse momentum k T . Directions of the lightcone variables 
fc ± are also indicated. In this figure we indicate, as black numbered points, momenta of 
three example emitted partons. Available phase space is limited from below by kf in = A, 



9 A sharp minimum/maximum rapidity is essential for avoiding mismatch between the parton distri- 
butions from two independent constrained MCs "operating" in the backward and forward hemispheres 
for the initial-state radiation. 
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Figure 4: Sudakov plane parametrized using two different sets of variables: (a) (t,v) and/or 
(77, lnk T ), (b) (i,v), where i = t — lnA , see Appendix. For simplicity x = 1 is set. Dashed 
red line marks position of the Landau pole. 

while in the direction k + the boundary is controlled by the total available q + , which is 
diminished by the factor a?j_i at z'-the step (setting for simplicity xq = 1). The minimum 
rapidity 77 = limits the emission phase space from the right hand side. The triangle 
and trapezoids show integration domains of the consecutive form-factors 
in eq. (1T4"|) , see also section 13.31 and Appendix IA.4I 

The above mapping of the evolution variables into four-momenta, with the minimum 
k T and built-in angular ordering, is identical to the one used in refs. [231 I2U [25], up to 
the following redefinition of the evolution time: tj = In Eh — r\i + In Xi-i = U + In Xi-%. For 
the redefined t we have 

i + In x > i n + In x n > i n -\ + lnx n _i > ... > £2 + In ^2 > ti + ln^i > t x , 

and consequently qi > with Zi = Xi/xi-i. 

3.2 Evolution kernels 

Once the phase-space parametrization is explained, we may define our choices for the 
evolution kernels, introduced so far only in the generic form in eq. (fT5|) . The kernels to 
be used in the CMC in this work will be of three kinds. The main difference between 
them is in the choice of the variable used as an argument of the coupling constant as- 
The appearance of the Landau pole in as will limit the choice of the IR cut-off in the 
multigluon phase space. Let us define first the gluonstrahlung kernel / = where / is 
flavour type, f = G,q, q, in all three cases. The strong coupling constant will be always 
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taken in the LL approximation 

= W- (26) 

s p lng-lnA 

Case (A): The standard DGLAP LL of ref. [8], which is used here as a reference case: 

Kf f <?(t t x,u) = ^ l -Pf){z) 6^ z>e = ^l l -Pfj{ X /u) 6 U - X>ue , (27) 

where e is infinitesimally small and z = x/u. 

Case (B): The argument in a$ is (l—z)e t = k T ju\ such a choice was already advocated 
in the early work of ref. [26]. For the IR cut-off we use A(t, u) = Awe - *, where A > Ao- It 
cannot be infinitesimally small, however, it becomes very small at large t. Using z = x/u 
we define 



>(0) 



ff \"i "7 - 1 f f V"V "1- z>\e~ t 

(21 



Case fCy 1 : The coupling constant 05 depends on the transverse momentum k T = 
(u — x)e t , while for an IR cut-off we choose A(t,u) = A(t) = Ae~*. The kernel reads: 

4 ( / C) (t,x,n) = a * {i \- X)e) 1 -Pj° f \x/u)9 u . x> _^, (29) 

For the (B) and (C) cases the choice of A for the IR cut-off must be such that we avoid 
the Landau pole in - the insertion of the 1 — z factor to the argument of the coupling 
constant resums higher order effects in the bremsstrahlung parts of the evolution [26] . In 
the non-diagonal, flavour-changing, elements of the kernel, there is no need for this kind 
of resummation. We can, therefore, use as(e t ). On the other hand, although quark gluon 
transition kernel elements have neither IR divergence nor a Landau pole, it makes sense 
to keep the restriction u — x < A(t, u), because the cut-out part of the phase space is 
(will be) already populated by the k T distribution of the primordial partorj^l. 

3.3 Virtual part of the kernel and form-factors 

For the evolution with any type of kernel the momentum sum rule 

= dt J dx xDf(t, x) — J du j^X^ J dx xJCff>(t,x,u)^j Dfi(t. 



u 



J du <|— uJC V ff(t, u) + J dx xJC 9 ff,(t, x, u)^j Df(t, 



(30) 



10 Without this restriction in the early evolution time quark gluon transitions would completely take 
over the bremsstrahlung. 
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is imposed, the same as in the reference DGLAP case. 

This sum rule determines unambiguously the virtual part of the kernel for all cases 
(A-C) 

*7/(*.«) =E J (31) 

It should be stressed, that in case (C) JC v (t, u) includes implicitly 9 u> A(t), as visualized in 
Fig. EJ^b). The following Sudakov form-factor results immediately: 

$ / (t 1 ,t |n)= / dt JC v ff (t,u) = / dt I —xlC e flf (t,x,u) 

J t f , J t Jo u 

, i (32) 
= / dt I — (u — y) ICfj-i f(t,u — y, u) = / dt J dz uzK, 6 f,f(t,uz,u), 

yv Jt Jo u j, J t Jo 

where z = x/u and y = u — x = (1 — z)u. The JCfi* is constructed using the IR- 
singular, non-singular and flavour-changing parts of the DGLAP (LL) kernel according 
to the following decomposition 



Qu-x>A(u)- (33) 



For the list of the coefficients Aff and the functions Ffrf(z) see Appendix of ref. |17j . 

We also need to define the generalized kernels beyond the case of bremsstrahlung, that 
is for the quark gluon transitions. One of the possible extensions, valid for all three cases 
X = A,B, C, reads 

xK?ffi(t,x,u) = S ff xK?ffi(t,x,u) + (1 - 5f*f) aS ^ ^ Fff(z)9 u _ x>A (x ){u) , (34) 

where as in the flavour changing elements has no z- or /c T -dependence and the IR cut-off 
is the same as in the bremsstrahlung case. The Sudakov form-factors resulting from 
the above kernels are split into three corresponding parts: 

$/(MoM = *o|«) + $ b f(h,t \u) + mt!,t \u). (35) 

In case (C) we get three genuinely w-dependent components 



& f (ti,t \u)= r dt i dz z)ue ) Aff e (l _ z)u>xe - t 

J to Jo 71 i — z 

2 

= Af f — g 2 {t + lntt,ti + lnu;t x ), 

Po 

rti rl n-y _ z ) ue t\ (36) 

$}(t 1 ,t |«)= / dt / dz-^ >- >-F ff {z) 9 {1 _ z)u>Xe - t , 

Jt JO 7r 

^(ti,t |«)= f 1 dt^->-J2 [ dzF n (z) e {1 _ z)u>Xe - t , 

J t 71 f'jtf'' 
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where U = ti — In A , i\ = t\ — ln A and function g 2 is defined in Appendix IA. 21 in terms of 
log functions. The integration domains for the consecutive form-factors of the above type 
are also shown in a pictorially way in Fig. m using the well known logarithmic Sudakov 
plane. 

In case (B) the u dependence disappears due to the fact that both — z)e t ) and 

@i- z >A(t) depend on u exclusively through z = x/u: 

*/(M„) = */(*!, *o|l), *<)) = $ 6 f (tl,t |l), mtlM = m^M^)- (37) 

In other words, setting u = 1 brings us from case (C) to the case (B). 

The reason behind the seemingly at hoc three-fold split of <&f(ti,to) is practical. In 
the MC the form-factor has to be calculated event-per-event. One-dimensional integration 
for each MC event is acceptable - it does not slow down MC generation noticeably. Here, 
the most singular part of is calculable analytically. In $^ we are able to integrate 
analytically over t and the integration over z is done numerically while in $j we can 
integrate analytically over z and the integration over t is has to be done numerically 
(see also Appendix IA.3I) . Altogether, we are thus able to avoid 2-dimensional numerical 
integration for each MC event (or use of look-up tables and interpolation for fast evaluation 
of the Sudakov form- factors for each MC event). 

Finally, the form-factors for the simplest DGLAP LL case read as follows: 

r(t )) A ff In-, 





2 

Wo 


(r(h 


*o) = 


2 

Wo 


(r(h 


h) = 


2 

Wo 





r(t )) J dz F ff (z) B x _ z>t , ( 38 ) 
r(to))Y] [ dz F f , f (z) 9^ z>e , 

tl-LtJO 



where r(t) = ln(t — In A ). 

At present, in the MC implementation, we use for cases (B) and (C), a slightly different 
form of the quark gluon changing kernels elements: 



'l-z>Ae- 



ir e ( B ')u \ x ir e ( B )u \ i ft x \ a s((l - z)e f ) 
Xitff {t, x,u) —Off xK,f)f[t, x, u) + (1 — Off) 1 ff\z 

xJC f ) f '(t,x,u) = 5ff xK,f)f'(t,x,u) + (1 - Off) Fff(z)9 y> xe- 



(39) 



that is, we use the same arguments of as for gluonstrahlung. We will refer to them 
as cases (B') and (C). This is done mainly to facilitate numerical comparisons with our 
Markovian MCs. One can easily go back from cases (B') and (C) to (B) and (C) with 
an extra (well behaving) MC weight, if needed. The corresponding form-factor <& C f(t\,to) 
gets properly redefined in cases (B') and (C), of course. 
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4 Constrained Monte Carlo for pure bremsstrahlung 



Let us discuss the case of pure gluonstrahlung first. We will focus our attention on the 
following integral, being part of eq. (fT4"]l 



G n ff (K B ;t b ,t a ,x,u) 



n rh 



[ / dU / dxi Xf f (ti 



° u x=x n ■ 

(40) 



It describes the emission of n gluons. The following "aliasing" of variables is used: t n+ \ = 
tb, t = t a and x = u. The integrand is well approximated by the product of the IR 
singularities in terms of variables y^ = Xi — 



n n _ n 1 n _ 



(41) 



Hence, switching from variables x$ to yi or Z{ is almost mandatory and the multigluon 
distribution with the 5-function constraining the total energy of emitted gluons takes the 
following symmetric form: 



®x=x r , 



(42) 

Constructing the MC program/algorithm for the multidimensional distribution featuring 
such a ^-function is a hard technical problem and it is the problem of constructing a 
Constrained Monte Carlo (CMC). 

It should be kept in mind, that it is possible, as shown in ref. [9], to generate the 
above distribution in x-space without such a 5- function, provided Gff is convoluted with 
the power-like function Xq . Such a solution was labeled as CMC class II, while the CMC 
of this paper was already referred to in ref. [9] as CMC class I. 

In ref. [8] the first CMC algorithm of the class I was found and tested for DGLAP 
kernel, that is for our case (A). This algorithm is based on the observation that for 
the product of steeply rising functions, proportional to the 5-function constraint 

is effectively resolved by a single (let's say) yk, while all other yi, can be considered as 
unconstrained. We shall extend the CMC class I solution to more complicated kernels, 
that is of our type (B) and (C). The main complication with respect to case (A) is due 
to a more complicated singular y-dependence (or ^-dependence) entering through the 
coupling constant as, and even more important, through form-factors. In addition, our 
new CMCs will not only generalize the solutions of ref. [S], but will be described in such 
a way that any future extension to other types of evolution kernels will be rather easy. In 
the following sub-sections we shall present the details of our generalized solutions. 
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4.1 Generic CMC class I 



As already indicated, we intend to introduce the formulation of the CMC algorithm which 
covers three types of kernels (A-C), and also that further extensions are possible and easy. 
At first, let us consider the following generic expression including the sum of constrained 
multidimensional integrals 



D(v) = V(v) 



+ ~\ / II K ( v i) dv i \ ■ (43) 

n=1 Uv i = i J ) 



The following properties will be assumed: (a) The function \&(i>), used to change inte- 
gration variables, must be monotonous; its derivative must remain non-negative, ^ = 
ty'(v) > 0, for all v > vq. (b) The positioning of the IR term0 5{v — vn) at v = vq is a con- 
vention. In practical MC realization it represents a no-emission eventl 12 !. Let us stress that 
in the w-space we are free to place the position vjr of the IR part of spectrum S(v — v IR ) 
anywhere outside the (v ,v x ) interval; we have opted for v IR = v (c) the variable v x 
controlling the overall normalization will be specified only later; it will be adjusted to get 
convenient normalization, (d) the true upper integration limit of Vi is below v and is in 
fact uniquely determined by the 5-function of the constraint. 

In the following examples, the variable V{ will be defined as Vi = In?/, = \iiXi — Xi-i or 
Vi = ln(l — Zi) = ln(l — Xi/xi-i), while the function ty(v) will be typically rather simple; 
= Vi or ^J(zi) = lnzi. 

Assuming K(v) > 0, we can define a mapping (and its inverse) which removes K(v) 
from the integrand: 



r< = R(vi) 



K(v')dv', Vi = V(n) = iT 1 ^). 



Our master formula transforms then as follows: 



r 00 1 r r R ^ n 1 1 

D[v) = m\v)e~ R{Vx) I 6 9 ( V)=9M + ~\ / II dTl <W)=£"=i m- 1 ^)) \ 

< 71=1 ' <-J0 i=1 J J 



(44) 



(45) 



The function f(r) = \I/(-R _1 (r)) is usually a very steeply growing function of r, hence 
the constraint is effectively resolved by a single rj ~ R(v), the biggest one. Let us 
exploit this fact in order to replace the complicated constraint with the simpler one^l 
S(R(v) — maxjTj). This is done in three steps. Step one: introduce a new auxiliary 
integration variable X countered by a 5-function: 



D(v) 



V=VQ 



X 



00 1 

Y- 

n=l 



^'(v)e- RM 



dX 



L Jo 



\dr t 



i=i 



5 *M=E"= 1 ^(«- 1 (^))^)= 



(46) 



11 This term is up to a constant (choice of integration variable) equivalent to S^,( v )—^( Vo y 

12 The "no-emission" MC event will have precise interpretation, independently of the choice of vq. 

13n 



The function max., rj is equal the biggest rj among j = 1,2, 



..,n. 
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Step two: change the variables = r[ — X: 
D(v) = e - RM 5 v=V0 + ^'(v)e~ R{Vx) 

-R(v) « 



dx / n^^>o 

n=l ' J L^ i=l 



5*(«)=E7= )-X))5R{v)-- 



(47) 



From now on we have to watch out for the condition = r\ — X > explicitly^- Step 
three: eliminate the old constraint by integrating over X 



00 i r rR(v) n 
D(v) = e- R ^S v=V0 + e~ R ^ £ - / \dr> 9 n>0 



" R(v)=ma,Xj r 



where Xo(r^, r^, r^) is found by means of solving numerically (iterative method) the 
original constraint for X. The Jacobian factor 3 is defined as follows 



(49) 



In the above factor the j-th term satisfying R(v ) = max,,- r'- dominates the sum in the 
denominator; consequently Z — R'{v) = K(v). This form, valid in the limit only, we keep 
explicitly in the integrand, while the remaining part of 2 is incorporated in a complicated 
but mild Monte Carlo weight: 



00 i r rR(v) n 

D(v) = e- R ^5 v=V0 + R'(v)e- R ^ ]T - / \dr> & 

n=l n -lJo i=1 J 



r' l= R(v) U!*, (50) 



where the MC weight is 



w 



# 



>0- 



(51) 



i=l 



Last but not least, we isolate cleanly a part of the integrand normalized properly to 
1 with the help of the integration variable £j = Ti/R(v) and the Poisson distribution 
P(n\X) = exp(-A)A n /n!: 



n „i 



D(v) = e~ R ^5 v=V0 + R'(v)e R ^- R ^ £ P(n - l\R(v)) [ / d& 



n=i 



i=l 



^maxjgj — 1 

n 



(52) 



The nice thing is that D(v), obtained by neglecting is known analytically 

D(v) = D(v)U =1 = e~ R ^6 v=V0 + 6 v>V0 R'(v)e R M- R M (53) 
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On the other hand, no problem with r[ < v, because we shall get X > 0. 
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and is normalized to 1: 




/ 

J c: 



1 



) 



D{v)dv = e 



+ 



d e R(v)-R(v x ) 



= e 



+ 1-e 



R{v x ) 



1. (54) 



cxp(—R(v x )) 



The above convenient unitary normalization is achieved by means of identifying v x with 
the upper limit of v, see also below for particular realizations. In the implementation of 
quark gluon transitions using FOAM [27] (see section \5§ one introduces the integration 
variable U = U(v) = e. R ^~ R ^ Vx ' G (0, 1) and the above equation transforms into 



The MC procedure of generating the variables v, n and (v 1 ,v 2 , ...,v n ) obeying the 
constraint is the following: 

• Generate v according to D(v) times whatever the other function of v in the Monte 
Carlo problem, for that purpose use the variable U = e R ( v )- R ( v ^) g (o, 1). 



• Otherwise U is translated into v (v > vq) and n > is generated according to 
P(n-l\R(v)). 

• Generate the variables £j G (0, = 1,2, ...,n, except one of them: ^- = 1, where 
j = 1, 2, ...n is chosen randomly with the uniform probability. 

• The variables £j are now translated into r[ and the transcendental equation defining 
the shift X is solved numerically. 

• Once X is known, then all ^(^(^((Ci))) are calculated. 

• The MC weight w# is evaluated; the check on > can be done earlier. 

The above algorithm generates weighted MC events (n; V\, v 2 , v n ) exactly according to 
the resumed series of the integrands defining D(v). 

4.2 Treatment of ^-ordering — generic case 

In case (A) of the DGLAP kernel, discussed in ref. [S], the t-ordering in the integrals 
of eq. (140]) can be easily traded for a 1/n! factor, while K(v) in the previous section in- 
cludes un-ordered integrals J dU over the entire available range for each t{. However, 
even in this simple case one has to be careful if one aims not only at the numeri- 
cal evaluation of the PDF, but also at the proper simulation of the entire MC events 
{n, (ti, xi), (t 2 , x 2 ), (t n , x n )}. The CMC class I algorithm of ref. [S] (DGLAP) is for- 
mulated in terms of Zi = Xi/xi-i, the point is that at the end of the MC algorithm one 
has to calculate Xi using a well defined ordering of Zi, which is exactly the same as in 




(55) 



• If U < e R ( Vx ^ then set v = vq and n = (no-emission event). 
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the sequence of the (originally) un-ordered tj. In practice, in the end of the MC gener- 
ation, one has to order U and Zi simultaneously before calculating Xi = XqZ\Z2---Zi. This 
is because the original integrand is symmetric with respect to interchange of the pairs of 
variables (tj, z^) <-» (t k , and not with respect to interchanging t{ t k or Zi <-» done 
independently. This property we will call the pairwise permutation symmetry. 

In our most general case (C), the integrand is not pairwise symmetric, mainly due to 
nontrivial x^-dependency in the form-factor, see eq. ( flOl) . In this case, the strategy is such 
that we introduce a simplified integrand, with the simplified form-factor and simplified 
kernel P (the IR-singular part of the kernel), such that the simplified integrand does 
feature the pairwise permutation symmetry. The above simplifications are immediately 
and exactly compensated by means of the MC weight u> p , which absorbs all possible 
pairwise non-symmetry. 

Let us translate what has been said above into rigorous algebra. We start from a 
more sophisticated variant of our generic multi-integral (I43p . covering all cases (A-C) and 
possibly other similar cases, but this time including explicitly ordered t-integrations: 



T^U 4. I \ <x,lt \ - ft d,t' f Vx (i „ TP(v' ,t')dv' J r 

D(t,t \v) =^'(v)e Jt o J "o(*') <5*( V)=9M + 



00 r n ^ pv 

+ Y1 [ dti dvi pfavi 

n=l L •/*»—! Jv (ti) 



U> P (t, V) 



(56) 



where t n = t. The bold-face variable t denotes the entire vector (ti,t2, ■■■,t n ) and similar 
convention is used for the definition of the vector v. 

In order to get rid of the t-ordering in the basic MC algorithm we proceed carefully 
step by step as follows: 

1. We introduce formally pairwise symmetrization, i.e. we sum up over all permuta- 
tions ^3 of the pairs of the variables (vi,ti), compensating by means of the 1/n! 
factor: 



r n rt 



n \^2 n , 



dt 



V0(tf) 



dvfntf.vf) 



F (tV 



2. The part of the integrand in the brackets is now perfectly pairwise symmetric and 
the permutation *P can be undone in this part: 



1 



HI 

i=l Jt ° 



dvi P(ti,Vi) 

vo(U) 



tf >tf, >...>* 



3. The sum over permutations J2<$ can be dropped out because only one permutation 
contributes at a given point t = (£j, £2, t n ). This particular permutation we denote 
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by ^t, obtaining: 

7-1/4. j I \ lT ,// \ -// dt' P{v',t')dv' J , 



oo r n 



+ J2 [ dU dvi P(ti,Vi) 

n=l L i=l •/"o(ti) 



(57) 



4. In the basic MC the weight u> F is temporarily neglected and U are generated un- 
ordered. The permutation <p t is then read from the ordering in t = (t\, t 2 , t n ). It 
is then used to construct the sequence of Xi out of Vi and to calculate w ¥ '. 

In order to finally bring eq. (157|) into the standardized form of eq. (1431) we interchange 
the order of integration over ij and 



t pv i>v rt( v i) 

dU / dviP(t h Vi) = / dvi / dtiP(t u Vi)— . ^ t ^ yviJ 

U-i Jvo{U) Jvo Jti-i(vi) Jvo JO 



dvi K(vi) / d(Ji. 



The additional change of variables 



a(t,v) 



r*™^ df p(f, 



K(v) 



(58) 



allows to map, for every value of V{, uniformly distributed <Ji into ij. For a particular 
realization see Appendix IA.3I Luckily, crj(tj) can be inverted analytically, i.e. ^((Tj,^) is 
available as an analytical formula, for all cases (A-C). 

The main purpose of the sub-section, was to obtain the following, ready for the MC 
implementation, representation of our standardized generic formula: 



D(t,t \v) = m'(v)e 



s:*K(v>)dv> 



i r n pt rl \ 

+ E T [ / dt * K ^ / da * 5 *W=E7- 1 ^ P (t^, v^) 
„=i n - \-i=i J u-i Jo J J 



i 



e-^)^ + ^ /2'(^)e^-^) 

oo p n „i „i 

xJ2P(n-l\R(v)) [ fa dat 



n=l 



i=l 



JO 



n 



(59) 



In the following subsection we shall describe three realizations of the above CMC class I 
schemes for three types of the kernels, (A-C). 
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4.3 CMC case (A), DGLAP 

Let us start with the particular realization of the above CMC for the easiest case of the 
DGLAP kernel, type (A). Such a CMC was exposed in detail already in ref. j8]. Here it 
serves as a reference case and a warm-up example. We recall the pure bremsstrahlung 
evolution operator of eq. (fl4l) in a form adopted for further manipulations: 



00 r n ft 



G ff (K B ;t bl t a ,x J u)=e-' s 'f^6 x=u + Y, U \ dU 6 

X e -*/(*i»*n) 



n=l i=l " la 



i=i 



(60) 



$x= 



where we have dropped the non-existing dependence on x$ in the form-factor $. Hence, 
one may exploit the relation Y17=o ^K^> U-i) = ^(^ t a ). After identification of terms and 
change of integration variables 



a s (e 



7i 



-ZiPff(Zi)9i- Zi>e , Zi = Xi/xi-i, 



a simplified expression is obtained: 

-G ff {K B - t b , t a , x, u) = e-*f^5 x=u + 
u 

00 r n „f b 
n=l L i=l "'^-i 



dt 5 - 



dz* ag 6 ZiPff(zi)9i- 



x/u 



7T 



■2i>£ 



(61) 



In the next step, the simplified kernel is introduced 

l-e l-e 



i=l ^ i=l ^ * 

n „ n „ 

= JJ / dvi as{e u )A f f = JJ / 

i=l, i=l„. 



(62) 



P(ti, Ui 



7T 



Vi = ln(l-Zi), w = lne. 



Once v is defined, we are ready to deduce our kernel K{y) and constraint function ty(v) 
from 



G ff (K B ;t b ,t a ,x,u) = e-^ itb ' ta) -5 Hx/u)=0 



00 r n 



X 

t b ln(l— x/u) 



+ e 



«e n * 



n=l L i=l 



du* P(*»,Ut) 



.r 



<5ln(x/u)=£™ ln(l-exp(^)) wf, 



ti-i 



I'd 



(63) 
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where 



n 

i=l 



a s (e u )zi{l - Zi)P f f{zij 



P(ti,Vi 



By means of comparison of the above expressions with eq. (1561) , see also eqs. 
can identify the components: 



v = \n{l-x/u), V(v) = ln(l - e v ), \V'(v) 



l-e v 



u — x 

X 



K(v) = T dt P(U, Vi) = A ff ^-(r(t b ) - T(t a )) . 

Jt a PO 

The expression for the basic form-factor can be identified also: 



(64) 



one 



(65) 



R(v) 



dv' K(v') 



dv' / dt'P(t',v') = A ff —(r(t b )~T(t a ))(v-v ). (66) 

Po 



The relation v = ln(l —x/u) is valid for v > Vq only. The variable v x represents the upper 
boundary oft' = ln(l— x/u), hence or fixed x and maximal u = 1 we obtain v x = m(l — x). 
With all the above elements at hand, we are able to complete the following standardized, 
accordingly to conventions of eq. fl59l) . formula 



-G ff (K B ;t b ,t a ,x,u) = \e- RM 5 v=V0 +6 v>V0 —^—R'\ 



v e 



R(v)-R(v x ) 



J2P(n-l\R(v)) ' d£i do { 

n=l \-i=l J ° Jo 



11 



(67) 



where 



w 



R{v x )-s f {t b ,ta) "Q ag(e f< )^(l ~ Zi)Pff(zi) 



(68) 



and the second weight is fully determined from eq. (jSTj) . supplemented with ^f(v) and 
K(v) of eq. (jB5jl. 

In the CMC, we usually integrate over u for fixed x, hence the following formula is 
relevant 

2 /• 

h<cxp(-R(v x ))\v=v + 9u>cxp{-R(v x )) 



J duG ff (K B ;t bl t a ,x,u) = dU (^j jfl. 

n „i 



n=l 



^max, £, =1 



11 



#(0 ^(t^v^ 



(69) 



where^l [/=[/( 



3 i?(«)-_R(^) 



e H(in(i- a: / u ))-fl(in(i- a: ))_ Due to re i a ti on R( Vq ) = 0, the 

point v = v o (IR boundary) translates into Uo = e~~ R ( Vx ' = e~ R ^ n ^~ x ^ . With all elements 
needed in eq. (|59|) at hand, we are ready to reconstruct from our generic formulation the 
complete CMC class I algorithm of ref . [8] , at least for pure bremsstrahlung. 



15 Another ingredient was the identity §M = (f^p 1 dR - L42 



du 
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4.4 CMC case (B), a 5 (e*(l - z)) 

This case is to some extent similar to the previous DGLAP case. We shall, therefore, 
concentrate on the differences. The starting point is now 



-G ff (K B ] t bj t a ,x } u) = e-*f^ t ^6 x=u + J2 HI dU / dx 

x e -*/(*i»,*»|i) 



oo r- n „t 



n=X L i=X 



TT— ^(ti.^x^Oe-^^W 

1=1 



(70) 



where one may again combine form-factors into a single one: ^" =0 = $>(tb,t a \l). 

Now, in terms of 



Xi%ff(ti, 



a s {e tl (l - Zi)) 



7T 



the simplified expression reads 

X -G Sf [K B - t b , t a , x, u) = e -*/&.*.|i) !^5 X=U + 

a s ((l - Zi)e li ) 



oo r n 



n=X L i=l 



+ [ dU dz 



ti-i 



x/u 



7T 



ZiPff{zi)e 



l—Zi>Xe 



Jx=u n"=i «j 



(71) 



Simplification of the kernel goes as follows: 



X-Xe-H 



II / dz * ^d 1 - ZiWziPffM - I] / _^L_ a 5 ((l - ^e**)^ 

i=l ^ t=l ^ ^ 

n „ n „ 

=n / ^(e^)^// = n / ]p(*i,«i), 



(72) 



where 



a s (e 



U+Vi 
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-Aff, Vi — ln(l — Zi) and Vo(t) — In A — t. 



(73) 



So far we have followed closely the DGLAP case, except of the more complicated IR cut- 
off and the extra factor in the argument of as- The choice of the variable v is the same, 
hence also the function ty(v) remains unchanged: 



G ff (K B ;t b ,t a ,x,u) = e-**&'*»>i{<J K i_ (j0)= o+ 

x I 



oo r n 



n=X L i=X 



( ^ln(l-e«)=E™ =1 ln(l-e , 'j) W X f) 



(74) 



U-i vo(U) 
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where 



P _ TT "5l« -jz^i - Zj)£ff{Zi) 
1 t\ ' 



While comparing the above expressions with eq. (J56J) we immediately identify the following 
components: 

P 7/ — T 

w = ln(l-x/n), V(v) = ln(l-e v ), \*'(v)\ = —— = ^-A (76) 

1 — e v x 

In principle, the evaluation of -fT(i>) requires a change of the integration order used in the 
form-factor integral, see also Appendix IA.31 

if, v v t(, ?j 

R(v;t b ,t a ) = J d£ J dv'P(t',v')= J dv' J dt'P(t',v') = J dv' K(v'). 

t a InA— t' fo(tf)) max(i a ,ln A— v') 'wo(ifj) 

(77) 

In practice, it is slightly easier to calculate R(v ; t b , t a ) with the t-integration external, and 
then obtain K(v ) by differentiation: 

h v 

2 f f 1 2 

R(v;t b ,t a ) = A ff — dt / dv'—— — — - = A ff — g 2 (t b + v,t a + v;t x ), (78) 

Po J J t + v — m Ao po 

t a InX-t' 

and 

2 

= d v R(v;t b ,t a ) = A ff —d v g 2 (t b + v,t a + v;t x ). (79) 

Po 

The functions q 2 and d v g 2 are given in appendices IA.2I and IA.31 The rest of the algebra 
is almost the same as in the case of DGLAP: 



^G ff (K B ] t b ,t a ,x,u) = \e~ J 



X 

n=l L i=l 



'V=VQ 


,0 1 1 

+ v>V0 X y>(v) 


</ 


doi 


^maxj £j=l 


Jo 


77, 



m -i L.-_wo Jo J ^ J 



The variable = ln(l — x) and the weights u>* and u> F are defined in the same way as 
discussed already in the case (A) of DGLAP. One has to remember only that enters 
into as{e tt+Vt ). The final integral form coincides with that for DGLAP, see eq. (IB"9"|) . The 
important differences are in the definitions of the components, in particular, the function 
R(v) is more complicated here. 
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4.5 CMC case (C), a s (k T ), CCFM 

For this case we are dealing with the most general implementation of the generic solution 
defined in eq. (JT] 



G ff (K B ;t b ,t a ,x,u) = e-^^5 x=u 



00 r n ft 



d = l 



+ E [ / <«* Wi / dx^%f f (t 



=1 ^ t a 



Xi-l 



e -E" =0 */fe^-il^-i)5 



where we are aliasing the following variables: t = t a and t n = tf,. The 

kernel is again more complicated than the one used in the two previous sub-sections 



ff \ ii^ii "^i—l) — ^■i^ff\^i)^Xi-\—Xi>Xe~ t ii %i — Xi/X 



7T 



i-l- 



With the help of the following transformation of the integration variables 



we obtain 

x 



1 Zj 



n 

, u Xj-i— Xi>Xe- f i 
CtX ^ X i 

Xi—i Xi 



dy% 



Ae" 



G ff (K B ;t b ,t a ,x,u) = e-*fto» t -M6 x=u + 



% 1 



+ e n / * / - ^^¥^.(1 - 



n=l L i=l 



■E5U 



(82) 



The simplification for the kernels and form-factor, to be compensated later on by the MC 
weight, reads as follows: 



-$f(tb,ta\ 



1 n 



i=l 



dyi 



a s (yie u )A ff 



(83) 



n „ n „ 

/(Wall-*) J / du . Q, s ( e M^) A// = e -*/(*6,ta|l-«) J / dViP(ti 

i_1 lnA-t, i_1 « (*i) 



P(U, Vi) = a s (e ti+Vi )A ff , Vi = Info), «o(<) = hi A - t, 
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where $/ is that of eq. ( 1361) . With the above definitions we obtain 

]G ff (K B ; t b , t a , x, u) = e-*t^ x -^h[x, u)^5 cxp{v)=0 + 



x 
-( 

u 



n=l 



j=l 



+ J2 dti dViP{U,Vi) 



exp(v)=J2" =1 exp(vj) 



W 



?4) 



U-i vo(U) 



where 
and 



v = \n(u-x), ^(v) = e v , ty\v) = e v = u - x 



p 



i=i 



55) 



(86) 



h(x, u) = e f 



# f(t b ,t a \l-x)-& f(t b ,t a \u) 



The factocj h(x,u) is compensating for the deliberate use of $f(t b ,t a \u) in the weight 
u> F , with the aim of ensuring u> p < 1. Since the simplified kernel P is the same as in the 
previous case (B), the standardized kernel and form-factor are also the same: 



R(v; t b , t a ) = A ff — Q 2 (t b + v, t a + v; t x ), K{v) = d v R(v; t b , t a ), 

Moreover, the upper phase space boundary is also the same v x = ln(l — x), hence 

R{v x ;t b ,t a ) = &f{t b ,t a \v x ). 
The standardized formula for the MC reads as follows 

X G ff (K B ;t b ,t a ,x,u) = h(x,u)\e~ R ^5 v=Mtb) + 6 v>Vo{tb) ^R' (v)e R ^- R ^ 

5„ 



57) 



u 



X 



J2P(n-l\R(v)) [ dfc / da 

n=l \-i=l J J 



The final integral, defining the distribution to be implemented in the CMC, takes the 
following form: 

J duGff(K B ;t b ,t a ,x,u) = J dU —h{x, / 6 , c/<cx P (— .rc^)) U=« (t b ) + du>ex P (-R(v x )) 



r Tl p\ p ^ 

VP(n -l\R(v)) I ft da 

n=l Li=l 7 ° J ° 



X 



— W W 



n 



(0 w*(t*v* 



(89) 



16 



It will be efficiently dealt with by the general purpose MC tool FOAM [27], see the next section. 
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where0 U = U(x,u) = e R ^~ R ^ = e R(Hu-x))-R(in(i-x)) _ The weight w # ig eva i uate d 
according to eq. (15 lj) . All important differences with the previous cases (A) and (B) are 
hidden in the definitions/constructions of the components of the above formula. The only 
explicit difference is in the presence of the factor h(x, u). 



4.6 Summary on CMC for bremsstrahlung 
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a s (e u )A ff 


Eq.([MD 
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InA - 1 


a s (e ti+Vi )A f f 


Eq.([75D 


Eqs. ([7911781) 
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u — X 


Xi—\ Xi 


e v 


InA-t 


a s (e^)A ff 


Eq.flM} 


Eq.(j87D 



Table 1: For three types of the evolution kernel %^ x \ X = A, B,C list of components in the 
generic eq. ( 1591) for CMC. Variable v x = ln(l — x) is always the same. 



Before we extend our CMC to the complete evolution with quark gluon transitions, 
let us summarize the case of pure bremsstrahlung. As we have seen, all three cases of the 
kernels (A-C) are compatible with the generic formula of eq. ( 1591) . provided we identify 
(define) properly all components there. These components are collected and compared in 
Tabled] for all three cases (A-C). The appearance, in case C, of the extra factor h(x,u) 
should be kept in mind. 

5 Complete CMC with quark gluon transitions 

Monte Carlo simulation/integration of variables related to quark gluon transitions is man- 
aged by the general purpose MC tool FOAM [27\. It has to be provided with the user- 
defined integrand, the so-called density function. Arguments of this function have to be 
inside a unit hyperrectangle (or a simplex). Starting from eq. (fTU|) we are going to reor- 
ganize explicitly its integration variables, see also the scheme in Fig. [TJ paying attention 

17 In this case R'(v)/$'(v) = dR/du. 
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to the integration order: 
i 

D f (t,x) = J dx G ff (K B ;t,t ;x,x ) D f (t ,x )+ 



E E 

n.=l /„_i,. ..,/i,/o 



n 1 1 1 

/ rft fc / dx n G ff {K B ;t,t 



X 



" n 1 1 

II J dx' k 7f kfk i (t k ,x k ,x' k ) J dx k ^i Gf k lfk l (K B ;t k ,t k _ 1 ,x' k ,x k _ 1 ) 



(90) 



where f n = f. Here, and in the following, we understand that the following integration 
order for all multiple integrations^ is used. 

J^J J ddi = J da n J da n ^\ • J da 2 J doi\. 

Let us now re-parametrize the above integral into a form better suited for integration by 
FOAM, that is in terms of 3n+l variables: Uq G (0, 1) and U k , ct k , /3 k £ (0, 1), k = 1, 2, ...n: 



D f (t,x)= [ dU (x,x ) H( x \x ,x) [ dG° W^D f (t , 

J JG{t b ,t a ) 



X ) + 



N 



+E E 

n.=l /„_i,. ..,/i,/o 



JJ / - t k -i)da k 
fc=i n 



dG (n) W? 



G(t b ,t n -i) 



X 



1 

n „ 

II / d & (l-X fc )^ /fc i (t fc ,X fc ,4(/? fc )) 
fc=l / 



x / dU k ^{x' k ,x k ^) H^ x \x k ^x' k ) / dG^" 1 ) Wti 



G(t fe ,t fc _i) 



■D/o(*o,x ) 



where for all three cases X = A, B,C, the multi-differentials <iG are used 19 ! 

/ duG ff (K B ;t b ,t a ,x,u) = [ dU H {x \x,u) [ dGW G 

Jx JO JG(t b ,t a ) 



(91) 



(92) 



18 This is the same order as for the operator products and in the time-ordered exponentials. 
19 See eqs. ^ and O. 
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they are defined in the following wa; 



/ dG^ = Ou<exp(-R(v x ))\v=v (t b ) + 

JG(t b ,t a ) 

oo p n „\ „i 

y2P(n-l\R(v)) n / dii / da, 

n=\ U=1 7 JO J 



U>exp(-R(v x )) 



<)- 



n 



(93) 



dU I dG= 1. 

G(t h ,t a ) 



As in case of dG^ k \ the weight 



(94) 



depend on the t range (in the above definition tb, t a is used). For each occurence of dG 
and W G in (EI]), different t range is used, as is seen from the context. We additionally 
mark the difference for the later use. 

For our three example kernels we have 



H {A \x,u) = H {B \x,u) = [-) 



it 



H {c) (x,u) = - h(x,u). 

X 



(95) 



The summation over the number of quark gluon transitions n is also mapped into one 
of the FOAM variables, after truncation to a finite N. In fact, the precision level of 
~ 10~ 4 is achieved for N = 4, see numerical results below. Following the prescription 
of ref. [8] the sum X]/ n _, /j / can be reduced just to a single term. Also, as in ref. [8], 
the additional mapping to the Tj = r(ti) variables is done in order to compensate partly 
for the t-dependence of as (U) in the ^ kernels. The purpose of this mapping is to boost 
slightly the efficiency of FOAM. 

FOAM generates all 3iV + 2 variables (including n) according to the integrand of 
eq. fl9TT) . omitting from it the following MC weight 



W 



IT". 

i=0 



G 



(96) 



Let us stress, that the calculation of this weight can be completed only after gluons 
are generated for all n + 1 bremsstrahlung segments first. Also, as in ref. [Hj, FOAM 
treats distributions as continuous in Ui- variables, ignoring 5-like structure in x' k — x^-i, 
corresponding to the no-emission case. This is very important and powerful method. The 
5-like no-emission terms are replaced by the integrals over finite intervals in Ui, exactly 
as in eq. 
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The index k in dG^ reminds us of the type of the parton type used implicitly in the distribution. 
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Figure 5: The distribution of rapidity and log of k T from CMC for 2E h = 1000 GeV, A = 1 GeV 
and the kernel type (C); pure bremsstrahlung case. 

With all the above formalism at hand, we can formulate our CMC algorithm similarly 
as in the general LL DGLAP case: 

• Generate super-level variables n, ti x\ and Xi with the help of the general-purpose 
MC tool FOAM according to eq. (|9"T|) . and neglecting W = J] W?. The parton 
types fi, i = 0,1, 2, 3..., n — 1 are determined from / = f n according to prescription 
of ref. 0. 

• In the above, the number of flavour transitions (G — > Q and Q — > G) is limited to 
n — 0, 1, 2, 3, 4, aiming at the precision of ~ 0.2%. 

• For each z-th pure gluon bremsstrahlung segment the emissions variables are gener- 
ated using the dedicated bremsstrahlung CMC of section HI according to eq. (l9~3l) . 
The weights Wf are calculated. 

• Generated MC events are weighted with W = TJ . They are optionally turned 
into weight=l events using the conventional rejection method. 
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• The four-momenta fcf and gf are reconstructed out of evolution variables and az- 
imuthal angled^!. 

The above algorithm is already implemented in the form of a program in C++ and tested 
using upgraded version of the Markovian MC of refs. [28] and [TTj . The numerical results 
are documented in the following section. 
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Figure 6: The gluon distribution from CMC and MMC for evolution with the kernel (C), 
A = 1 GeV and e imax = 2E h = 1 TeV. Contributions from fixed number of the quark gluon 
transitions n — 0, 1, 2, 3, 4 are also shown. The ratios CMC /MMC in the lower plot are given 
separately for the total result and for the number of quark gluon transition n = 0, 1,2. The 
MC statistics is 4.5 • 10 9 weighted events for CMC and 10 10 for MMC. 



6 Numerical tests 

Most of the numerical tests proving that the new CMCs of this paper work correctly were 
done by means of comparison with the updated version of the Markovian MC program 

21 Azimuthal angles are generated uniformly. 
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log (x) 



Figure 7: The Quark distribution from CMC and MMC for evolution with the kernel (C), 
A = 1 GeV and e* max = 2Eh = 1 TeV. Contributions from fixed number of the quark gluon 
transitions n = 0, 1, 2, 3, 4 are also shown. The ratios CMC/MMC in the lower plot are given 
separately for the total result and for the number of quark gluon transition n = 0, 1, 2. The 
MC statistics is 10 10 weighted events for both CMC and MMC. 



MMC [15], which is a descendant of that described in refs. [171 [28]. The version of MMC 
used below implements exactly the same type evolution with exactly the same kernels and 
boundary conditions. We can therefore expect numerical agreement of CMC and MMC to 
within statistical MC error, which will be of the order of 10~ 3 , for about 10 10 MC events 
(and for about 10 2 different values of a; in a single MC run). 

We start, however, with the simple tests in which we verify the correctness of the 
mapping of the evolution variables into four-momenta. In Fig. [5] the distribution of 
rapidity and kj of the emitted gluons is shown. Sharp cut-offs corresponding to the 
minimum rapidity rj ^ (maximum f) and minimum k T ^ 1 GeV are clearly visible in 
the plot. Note that this plot shows the same triangular area of the logarithmic Sudakov 
plane as in Fig. BJa), but now populated with the MC events. 
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I Multiplicity 





Figure 9: The distributions of n and Xi and the emission multiplicity from the CMC and 
MMC programs. The results are from the MC run for quark in proton. 



6.1 Testing CMC versus MMC 

We examine results of the evolution from the initial energy scale A = 1 GeV up to final 
energy scale of Eh = 1 TeV, using exactly the same initial quark and gluon distributions 
in a proton at the Qq = e*° = A = 1 GeV scale as in ref. [8] (they were also used in 
ref. [2H]). 

In Fig. [6] the distributions of gluons at Eh = 1 TeV are shown, while in Fig. [7] the 
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corresponding results for quarks, Q = q + q, are exposed, also at the same high energy 
scale Eh. = 1 TeV. Numerical results are provided for the evolution type (C) (see section 
13.21) . that is for as{k T ) and the evolution time being identical with the rapidity of the 
emitted particle. In these two figures we compare the gluon and quark distributions 
obtained from the current CMC program and from the updated MMC program [T5] . 
The principal numerical results from both programs, marked as "total" in the upper 
plot of both figures, are indistinguishable. We therefore plot their ratio in the lower 
plot of both figures. The parton distributions from the CMC and MMC programs agree 
perfectly, within the statistical errors of ~ 10 -3 in the entire range of x. In the same plots 
we include individual contributions from a fixed number of the quark gluon transitions 
n = 0, 1, 2, 3, 4, and their ratios. Again very good agreement between CMC and MMC 
results is seerO 




It should be stressed that the MMC program has been tested numerically at the same 
~ precision level for all types of the kernels (A-C) by comparison with the semi- 
analytical (non-MC) program APCheb [29]. In addition, for the DGLAP case, MMC was 
also successfully compared with another non-MC program QCDNuml6 [30]. These dedicated 
tests of MMC using non-MC programs will be published separately in the forthcoming 
paper [15] . 

Additionally, in Figs. [8] and [9j we show comparisons of selected (semi-) exclusive dis- 
tributions from CMC and MMC, i.e. we have chosen for the test the distributions of the 
consecutive evolution time variables Tj(ij), the energy fractions Xj and the total parton 
emission multiplicity from both CMC and MMC. The above distributions resulting from 
separate CMC and MMC runs are interposed - no visible difference is seen within the 
resolution of the plots. The above plots are for a final parton being a gluon in Fig. [8] and 
a quark in Fig. [9j In these plots one is testing a nontrivial aspect of the CMC algorithm 
related to removing and restoring the t-ordering described in section 14.21 In fact, any 
departure from the pairwise ordering procedure described in section 14.21 would ruin the 
agreement of CMC and MMC seen in Figs. [S] and M The weight distribution from CMC 
is also shown in the lower part of these figures. 

All the above tests were done for the most interesting case of CMC with the kernel 
(C). What about numerical tests for cases (A) and (B)? In ref.[S] numerical agreement 
within ~ 10~ 3 statistical error between CMC and MMC for the LL DGLAP case (A) 
was already documented. While preparing this work, we have compared CMC and MMC 
for case (B), as(e _< A), obtaining equally good numerical agreement. The corresponding 
plotS look quite similar to these in Fig. [6] and Fig. 

22 Certain disagreements for n = 4, which affect total result at the 10~ 4 level can be traced back to 
well understood inefficiencies of FOAM at higher dimensions. 

23 They were shown in the contribution by S. Jadach to the Ustroh Conference, September 2005, see 
http : / /home . cern . ch/ j adach/public/ustron05 . pdf 
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7 Discussions and summary 



We have generalized the constrained Monte Carlo algorithm of ref. [8] from DGLAP to two 
other more complicated types of the evolution equations, one of them fully compatibly 
with the CCFM evolution equation [12]. The above extentions of the CMC algorithm are 
implemented in the computer program CMC, tested by comparisons it with the Markovian 
MC, which also solves exactly the same evolution equations. 

Let us point out that the most complicated version (C) of the evolution equation 
elaborated in this work follows closely the CCFM model as formulated in [25] (except for 
the temporarily omitted non-Sudakov form- factors) . Its Markovian MC implementation 
SMALLX was worked out in refs. [231 GI]> an d later on exploited in the construction of the 
CASCADE MC of ref. |31j . which is based on the backward evolution algorithm [7]. Our 
CMC program was tested against our own Markovian MC, see previous section. However, 
it would be interesting to compare it also with the above SMALLX and CASCADE programs. 
We hope this will be done in the near future. 

Since the CMC program of this work will be used as a building block in the MC event 
generator for the Drell-Yan type processes and deep inelastic lepton-hadron scattering, 
the explicit mapping of the evolution variables into four-momenta of the emitted particles 
is also defined, implemented and tested. The evolution time is mapped into the rapidity 
variable (angular ordering). A sharp cut-off is imposed on the k T and rapidity of the 
emitted particle. The sharp cut-off in the rapidity will be useful when combining two 
CMCs for two colliding initial state hadrons, with neither gaps nor overlaps in the emission 
phase space. 

The CMC algorithm was worked out in detail and tested for three types of the evolution 
kernels and the phase-space limits. It was purposely defined/described in such a way that 
it can be easily generalized to other types of the evolution kernels. In the following works 
it will be used as a component in the MC modelling of the initial state parton shower for 
showering of a single hadron in a more complete MC project for LHC. 
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24 Wc did not include into the discussion the so called non-Sudakov form- factor for CCFM, case (C). 
However, it is already included in the code. 
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Appendix 

A Auxiliary functions in form-factors 
A.l Triangle function g 




Figure 10: The integration area in the definition of the functions: (a) g(t b ;t\) and (b) 
g(t b + v x ;t x ). 

The simplest component function in the Sudakov formactor in cases (B) and (C) of 
^-dependent as reads as follows 

g(t b ;i x ) = j'dtj dv^t^ = 6t b>h {i b [lni b -lni x -l}+ix}, (97) 

where A defines the IR cut-off. Here and in the following, we follow certain notation rules 
allowing us to write the above and similar functions in a compact way: 

1. For all variables like t, t b and t\ bar over them means i = t — InAo, i b = t b — InAo, 
i\ = t\ — hiA , etc., where A is that in eq. ( 1261) . i.e. the position of the Landau 
pole. 

2. Occasionally we shall omit the explicit dependence on t\ = In A, that is we always 
understand g(t b ) = g(i b ;ix). 

3. We always understand that j^dx = when b < a. 

The area of the integration in eq. ( 1971) is the triangle, depicted if Fig. [10] (a). 
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The following similar integral, with the same integrand, but slightly different triangular 
integration area, depicted in Fig. [TO] (b), can be expressed using the same function g: 




Figure 11: The integration area in the functions: (a) g2(t a , h] tx) and (b) Q2(t a +v x , h+v x ; t\). 
The other basic IR-singular component function in the Sudakov form-factor reads 

ftj(ta, hi h) = [" di [ dv 9 -^^ = e h> t a {g(t b ] h) - g(t a ; h)} (99) 

Jt a J t + V 

The corresponding integration area is the trapezoid depicted in Fig. [TTJ (a). The above 
calculation result is obvious if one notices that the trapezoid is the difference of two 
overlapping triangles. Similarly, the following integral with the trapezoid integration area 
of Fig. [11] (a) is again expressed in terms of the functions already defined above 

Q2 (v x ,ta,t b ;tx) = / dt dv-j v - 



t + v (100) 

= Qt b >taie(ta + V x ] t X ) - Q(tb + V x ] tx)} = 02 (*a + V x , t b + V x ] t\). 

Let us stress that keeping 9t b> t x in the basic function g of eq. ( 1971) is essential for validity 
of the evaluation of all the following related functions and integrals. 
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IF 1 ]—— — 

Figure 12: The integration area in the function g l 2 (t a ,h,v x ;t\). 
A. 3 Non-singular functions and mapping 

Let us evaluate now the following integral, which is similar to Q2, except that we insert 
into the integrand the additional function F(v,t): 



e [ P(v x ,t a ,t b ;t x ) = / dt dv^±F(v,t). (101) 



6, 



Of course, for F — 1 we recover g 2 . However, even in this case the following evaluation of 
q\ f ^ will be interesting, because we shall swap the integration order and introduce variable 
mapping, exactly the same as we have used for finding out K(v) and eliminating K(v) 
through the additional mapping r{v). Such a swapping integration order is also used 
for evaluation of the non-IR part of the bremsstrahlung integral and flavour changing 
contributions, where we have F(v,t) = F(v), hence the inner integration over i can be 
done analytically and the outer one over v is done numerically. 
In a general case, after swapping the integration order 



Q 2 F] (v x ,i a ,i b ;t x ) = J dv{6 Vx<h . la / dt r f- + 9 Vx> t x -t a I dt r f- }F(v, t) (102) 

t\—tb 

the integration is split into two parts, first for the triangle and second for the rectangle 
integration area, see Fig. [12] for illustration. The above change of the integration order 
was also exploited by the authors of HERWIG MC [5j 132^1 . The internal integral can be 




3 We acknowledge private communication from Mike Seymour for on that, see also Chapter 5 in 



I http:/ /hepwww. rl.ac.uk/ theory /seymour/ thesis/ 
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transformed using the identities: 



o 



e'(t b ) = dt b e(t b ; t x ) = J dv = 9t b> i x (\nt b - lnt A ), 



t\ —th 



Qi (v x ,t a ,t b ) = d Vx Q [ 2 ] (v x ,t a ,t b ;t x ) 

= 8v x >t x -i a W(t b + v x ) - e'(t a + v x )\ + 6 Vx< i x -t a e'(i b + v x ), 

to the following convenient form 



i 



(103) 



8 [ 2 ] (v x , t a , h\ t x ) = I dv gf ] (v, t a , t b ) I d ( | n (^ + v ^ J ^ 

tx-tt ^ (104) 

= e2 ] (v x ,t a ,t b ) J d£(v) J da(v,t) F(v,t), 



where 

a \ S ] ( v Ja, t b ) i t\ ln(t + u) 

^^ = ^U7 — — 7' a M) = /[ii, - ( 105 ) 
Q 2 {v x ,t a ,tb) Q-2 {v,t a ,t b ) 

The above mapping is used in the bremsstrahlung CMC in cases (B) and (C). In this 
context one also needs to perform the inverse mapping t>(£), which requires numerical 
inversion of q\ 1 \v, tai tb) as a function of v. Once v is known, the second inverse mapping 
t(cr, £) is easily implemented, as it can be formulated in a fully analytical form. As already 
indicated, a very similar variant of the above integration order is used in the evaluation 
of the non-IR and flavour- changing parts of the Sudakov form-factor. 



A. 4 Integration area in the plane of In A; and rapidity ij 

Finally, let us relate the phase space v and t used in the calculation of the form-factor 
above with the Sudakov plane in lnk T and the rapidity rj. This is done in the pictorial 
way in Fig. HI where we depict a situation with three emission, underlining the trapezoid 
integration domain for the Sudakov function $(^3,^21^2)- This is show on the right hand 
side of the figure, as the trapezoid marked by abed, in the plane of t and v. On the left hand 
side of this figure, the corresponding trapezoid is seen in the plane of In k T and rapidity r\. 
Let us remind the reader that in the case (C) we define Vi = In?/, = \n(xi — Xi-i) and the 
rapidities t]i are related to the evolution times U with the simple linear transformation of 
eq. (J2ID in Sec. I3TH 
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